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1.  Introduction,  Summary,  and  Recommendations 


The  objective  of  this  project  was  to  explore  ways  to  process  magnetometer  or 
electromagnetic  induction  (EMI)  sensor  data  so  that  UXO  discrimination  processing  UXO  can  be 
enhanced.  Suppression  of  clutter  was  a  principal  focus.  To  this  end,  magnetic  field  data  obtained 
above  and  near  the  ground  surface  was  regarded  as  something  of  a  boundary  condition.  As  such, 
it  should  allow  inference  of  magnetic  fields  or  associated  source  distributions  either  upwards  or 
downwards.  The  former  is  done  to  reduce  noise  and  clutter  signals,  the  latter  to  observe 
concentration  of  field  or  source  values  around  the  location  of  a  subsurface  object  of  interest.  This 
was  done  here  solely  on  the  basis  of  the  governing  physics,  without  application  of  any  particular 
kind  of  target  modeling  and  its  attendant  optimization  searches  for  position,  orientation,  and 
target  properties.  Classical  geophysical  processing  of  this  sort,  using  the  data  over  a  measurement 
surface  from  a  single,  widespread  response  field,  constitutes  “continuation”  of  that  field.  To  our 
knowledge,  this  study  is  the  first  attempt  to  apply  the  concepts  to  transient  electromagnetic 
induction  surveying,  particularly  for  the  kind  of  high-resolution,  shallow  sensing  needed  for  UXO 
discrimination.  Our  application  in  EMI  is  distinguished  by  the  fact  that  oblique  (non- vertical) 
magnetic  field  components  must  be  used,  and  the  associated  source  distributions  computed  must 
be  displaced  from  the  plane  of  measurement. 

For  true  continuation,  it  is  essential  that  the  measurements  all  pertain  to  a  single,  though 
spatially  varying  field.  In  magnetometry  surveying,  the  excitation  (primary)  field  is  the  earth’s 
magnetic  field.  For  all  practical  purposes,  this  static  field  is  uniform  over  the  survey  area  around  a 
target.  All  measured  responses  therefore  constitute  static  responses  to  the  same  excitation.  In 
particular,  all  data  obtained  around  a  buried  object  simply  constitute  samples  at  different  locations 
of  the  same  response  field  emanating  from  the  target,  as  that  field  varies  in  magnitude  and  vector 
orientation  from  place  to  place.  The  time  domain  (TD)  EMI  system  used  in  this  project  is  similar 
in  its  nature,  despite  the  transience  in  the  responses.  Here  the  transmitted  primary  field,  while 
different  from  the  earth’s  field,  is  also  unchanging  in  magnitude  and  orientation  over  the  region 
around  the  target.  While  surveyors  rove  about  with  local  receivers,  the  excitation  of  the  target  is 
always  the  same;  its  internal  response  is  always  the  same;  and  it  is  always  producing  the  same, 
spatially  varying  exterior  distribution  of  magnetic  field  in  response,  which  is  observed  from 
different  vantage  points.  As  the  TD  EMI  data  obtained  here  constitute  a  distributed  sampling  of 
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the  same,  single  response  field  it  can  likewise  serve  as  a  consistent  boundary  condition,  or  the 
equivalent,  for  that  field.  The  receivers  used  here  are  magnetometers  of  much  the  same  kind  as 
are  used  in  static  magnetometry  surveys.  At  each  point  in  time  they  simply  catch  a  snapshot  of  the 
transient  magnetic  field.  All  results  developed  here  for  TD  EMI  thus  apply  directly  to  treatment 
of  magnetometry  data. 

Note  that  the  nature  of  the  data  described  above  is  fundamentally  different  from  that 
obtained  with  most  of  the  popular  man-  or  cart-portable  EMI  sensors  used  in  UXO  surveying  and 
related  studies  (EM-61,  EM-63,  GEM  sensors...).  Those  sensors  are  “mono-static”  in  the  sense 
that  local  transmitters  and  receivers  are  approximately  co-located.  As  a  sensor  of  this  sort 
proceeds  through  a  survey,  the  spatial  relation  between  the  transmitter  and  target  constantly 
changes.  The  primary  field  impinging  on  the  target  constantly  changes.  Data  obtained  at  each 
location  is  a  sample  only  of  the  field  in  response  to  the  excitation  transmitted  from  that  location. 
Therefore  a  spatial  array  of  such  data  cannot  constitute  a  boundary  condition,  or  the  equivalent,  in 
the  same  sense  as  used  above.  One  might  say  that  data  at  different  points  constitute  samples  of 
different  boundary  conditions,  i.e.  from  a  different  fields.  Consequently  the  data  cannot  be  used 
for  continuation,  at  least  in  the  classical  sense.  We  return  to  the  topic  of  upward  or  downward 
projection  of  data  from  mono-static  sensors  in  our  recommendations  below.  Meanwhile,  all  work 
in  this  project  pertains  to  more  classical  continuation  of  single  fields. 

To  our  knowledge,  in  connection  with  surveying  applicable  to  UXO  discrimination  only 
some  short  distance  upward  continuation  of  static  magnetometry  data  has  been  applied  heretofore 
[1].  The  classic  1965  text  by  Grant  and  West  [2]  still  provides  an  excellent  treatment  of  the 
fundamental  relations  needed,  as  well  as  some  vital  computational  maneuvers,  if  only  discussed 
with  reference  to  static  magnetic  and  gravitational  fields.  In  those  areas  of  application,  the 
possible  utility  of  upward  and,  occasionally,  downward  continuation  are  recognized  in  well 
known  contemporary  geophysics  texts,  e.g.  [3,  4] 

Clutter  is  the  bane  of  EMI  surveying  for  subsurface  UXOs  under  realistic  circumstances. 
Isolating  the  contributions  of  different  signal  sources  and  suppressing  clutter  responses  would 
take  us  a  long  way  towards  reliable  subsurface  UXO  discrimination.  The  work  performed  here 
successfully  suppressed  clutter  in  field  data,  even  when  the  local  clutter  signal  was  two  to  three 
times  that  of  the  broader  underlying  UXO  response  and  was  situated  directly  on  and  within  it. 

The  success  of  the  approach  is  directly  tied  to  the  fact  that  it  relies  on  the  governing  physics  [5, 
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6],  and  not  on  signal  processing  stratagems  divorced  from  that  physics  or  idealizing  it  unduly. 
Beyond  use  of  the  applicable  physics  in  the  numerical  formulations  and  computations,  that 
physics  is  also  at  the  root  of  the  phenomenology  exploited  for  advantage  in  upward  continuation. 
The  relevant  phenomenology  is  explained  briefly  as  follows. 

Consider  the  response  from  a  discrete  object  that  is  subject  always  to  the  same  primary 
field,  when  the  receiver  is  at  least  a  characteristic  length  from  the  object.  The  magnitude  of  the 
signal  is  (approximately)  inversely  proportional  to  the  cube  of  the  distance  d  between  receiver 
and  object. 


b 

~d^ 


(1) 


where  the  constant  of  proportionality  a  depends  partly  on  shape  and  object  composition  but  most 
significantly  on  the  size  of  the  object.  For  example,  in  the  case  of  a  sphere,  b  is  proportional  to  the 
object’s  volume.  Thus  a  sphere  with  a  diameter  that  is  substantially  larger  than  that  of  another 
will  produce  a  larger  signal  in  proportion  to  the  cube  of  the  ratio  of  their  diameters,  if  other  things 
are  equal.  We  assume  here  that,  from  the  point  of  view  of  UXO  clean  up,  to  qualify  as  a  clutter 
item  an  object  must  be  significantly  smaller  than  a  target  of  interest,  i.e.  UXO.  So,  why  would 
one  worry  about  an  inherently  smaller  signal?  The  answer  is  that  “other  things  are  [not!]  equal.” 
By  (1),  a  small  metallic  item  on  the  ground  surface  that  is,  say,  ten  times  closer  to  a  sensor  than  a 
buried  UXO  gains  in  relative  signal  magnitude  by  a  factor  of  10^,  relative  to  what  would  be  seen 
if  they  were  the  same  distance  away.  This  can  often  more  than  make  up  for  the  inherent  weakness 
of  the  clutter’s  response,  raising  it  to  problematical  prominence.  Overall,  clutter  will  only  be 
important  if  it  is  near  the  surface. 

Upward  continuation  can  be  used  to  take  advantage  of  this  kind  of  power  law  signal 
dependence  to  suppress  the  relative  importance  of  near-surface  clutter.  From  the  discussion 
above,  if  the  magnitudes  of  two  signals  from  a  given  object  are  compared,  when  obtained  from 
different  distances  d\  and  <^2,  the  ratio  of  signal  strengths  will  be 


yd^j 


(2) 
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Thus  raising  the  observation  point  relative  to  a  near-surface  piece  of  clutter  will  engender  a  much 
larger  proportional  diminution  of  signal  strength  than  would  be  the  case  for  a  deeper  target  of 
interest.  See  Figure  1.  As  the  sensor  is  raised,  the  quantity  (dci/dc2)^  diminishes  the  clutter’s 
response  by  a  much  greater  amount  than  (dti/dt2)^  reduces  the  UXO’s. 


Sensor 


Figure  1.  Two  elevations  of  a  sensor  (receiver)  relative  to  a  combination  of  near-surface 
clutter  item  and  deeper  UXO. 


This  effect  was  explored  in  some  detail  by  the  PI  in  another  SERDP  project  (MM- 1282) 
and  was  reported  in  its  Final  Technical  Report  as  well  as  reference  [7].  As  noted  in  those 
references,  the  principal  difficulty  in  trying  to  take  advantage  of  this  phenomenon  is  that  it  may 
be  impractical  to  raise  the  sensor.  It  may  be  physically  quite  difficult  to  do  so,  and  perhaps  more 
fundamentally,  all  the  signals  involved  may  become  too  weak  to  be  useful.  Basically,  in  trying  to 
shift  the  relative  balance  of  target  and  clutter  signals  one  may  produce  overall  SNRs  that  are 
intolerable.  The  motivation  for  pursuing  upward  continuation  (UC)  comes  specifically  from  the 
recognition  that  it  will  not  be  troubled  by  these  problems.  We  do  not  have  to  worry  about 
generating  worse  relative  noise  levels  at  greater  distances  if  we  do  not  actually  have  to  perform 
the  measurements.  More  specifically,  UC  will  actually  tend  to  increase  the  SNR  relative  to 
random  noise  and  error.  Just  as  UC  will  reduce  the  kind  of  signal  perturbation  caused  by  clutter 
such  as  that  indicated  schematically  in  Figure  1,  it  will  also  reduce  irregularities  caused  by  other 
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perturbations  such  as  uneven  ground  effects,  positional  errors,  etc.  Both  simulations  (Section  3.1) 
and  field  tests  (Section  3.2)  demonstrate  this,  while  Section  5  explores  some  of  the  basis  for  it. 

Downward  continuation  (DCN)  has  various  meanings  in  the  geophysics  literature, 
sometimes  phasing  into  broader  inferences  of  subsurface  structures  and  phenomena  than  just 
computation  of  signals  at  lower  (subsurface)  elevations.  In  any  case,  its  purpose  is  generally  to 
isolate  subsurface  sources  for  signals  of  interest.  Here,  UC  was  accomplished  by  solving  for  a 
distribution  of  equivalent  magnetic  sources  over  some  plane  at  an  arbitrary  distance  below  that  on 
which  data  were  obtained.  These  sources  can  then  predict  the  signal  that  would  be  obtained  at  any 
location  above  the  original  data  plane.  This  is  a  much  simpler  computational  approach  than  those 
in  noise-amplifying  differential  equation  UC  approaches.  Our  UC  system  goes  through  the  same 
computational  motions  as  would  be  required  to  solve  for  sources  over  deeper  and  deeper  planes, 
as  in  a  variety  of  DCN.  The  same  overall  algorithms  were  used  here  for  both  UC  and  (what  we 
are  calling)  downwards  continuation. 

Solution  for  sources  over  a  plane  displaced  from  the  data  plane  brings  ill-conditioning 
issues  to  the  fore,  for  reasons  explained  in  Section  5.  The  deeper  the  source  plane,  the  worse  the 
ill-conditioning,  until  a  very  modest  amount  of  noise  or  clutter  in  data  causes  completely 
unacceptable  distortions  in  the  computations.  We  show  here  that  this  ill-conditioning  is 
completely  controllable  by  means  of  the  spectral  truncation  methods.  Two  methods  are 
investigated,  namely,  1)  Fourier  transformation  of  both  data  and  source  distributions  into  the 
spatial  (not  electromagnetic)  frequency  domain,  with  truncation  of  higher  and/or  weaker 
frequency  content;  and  2)  transformation  into  an  eigen-system  framework  and  truncation  of 
equations  with  weak  eigenvalues.  The  gist  of  each  of  these  approaches  is  to  translate  the  problem 
into  a  basis  in  which  each  degree  of  freedom  (coefficient)  has  a  distinctly  independent  influence 
on  the  system,  and  weak  participants  are  readily  recognized  and  discarded.  The  results  suggest 
that,  while  both  methods  stabilize  the  solution  impressively  in  the  face  of  enormous  condition 
numbers,  neither  is  likely  to  be  very  useful  for  focusing  on  subsurface  source  locations.  The 
spectral  truncation  necessary  for  stabilization  limits  the  resolution  of  the  computations.  At  the 
same  time,  the  frequency  truncation  (FT)  approach  appears  superior  in  consistency  and  reliability, 
which  relates  to  the  nature  of  the  integral  operators  used.  Examples  suggest  that  FT  may  also  be 
able  to  suppress  amplification  of  equivalent  sources  associated  with  clutter  items  more  shallow 
than  the  plane  of  equivalent  sources.  This  strong  clutter  amplification  effect,  which  compounds 
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but  is  different  from  noise  amplification  by  ill-conditioning,  is  in  itself  enough  to  recommend 
caution,  if  not  complete  distaste  for  DCN. 

Conclusions  together  with  recommendations  for  follow-on  work  are  as  follows: 

1 .  Downward  continuation  (DCN)  appears  to  have  limited  appeal,  even  if  the  daunting 
problems  of  ill-conditioning  and  clutter  amplification  are  dealt  with  successfully.  Spatial 
frequency  truncation  (FT)  appears  capable  of  stabilizing  DCN  calculations  while 
suppressing  clutter  amplification.  FT  warrants  further  investigation,  particularly 
because  some  of  the  same  problems  appear  in  the  generally  more  appealing  upward 
continuation  (UC)  computations. 

2.  Upward  continuation  was  eminently  successful  at  clutter  and  noise  suppression  in 
all  simulation  and  field  examples  considered.  While  ill-conditioning  issues  can  arise 
even  in  UC  calculations,  they  appear  to  be  readily  controlled  by  spectral  truncation  of  the 
data  and  of  the  inferred  source  distributions.  Data  were  continued  easily  up  to  elevations 
where  signals  were  dominated  by  a  deeper  UXO,  even  when  very  strong  signals  from 
near-surface  clutter  items  were  present  in  the  original  data.  This  now  needs  to  be  tested 
in  discrimination  studies  that  exploit  such  continued  data.  Note  that  the  continuations 
entail  no  idealization  or  non-physical  simplifications  of  responses  predicted  for  higher 
observation  points. 

3.  In  conjunction  with  item  2,  sensitivities,  minimum  data  and  resolution  requirements 
need  to  be  determined  for  UC.  Here  only  very  modest  resolutions  were  employed,  with 
no  apparent  ill-effect.  Only  ~  10  cm  point  spacing  was  used  in  the  processed  data  arrays, 
with  underlying  cross-track  resolution  really  more  like  ~  20  cm.  Because  UC  signals 
were  generated  computationally  using  integral  incorporation  of  equivalent  sources,  only 
the  same  undemanding  resolution  was  required  of  the  source  distributions.  This  is  much 
simpler  and  requires  much  less  resolution  that  differential-equation-based  calculation  of 
upper  or  lower  fields  from  measured  boundary  conditions.  That  said,  surely  some 
resolution  constraints  will  apply.  In  the  examples  considered  with  a  sharp  clutter  spike 
directly  on  a  broader  target  signal  pattern,  the  success  of  the  method  obviously  depended 
on  the  system’s  ability  to  discern  that  the  clutter  produced  a  sharp  divergence  in  signal 
value,  i.e.  a  rapid  change  in  value  over  space.  Beyond  this,  fundamental  uniqueness  and 
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invertability  questions  regarding  equivalent  source  distributions  should  be  clarified 
beyond  what  was  undertaken  in  this  brief  study. 

4.  The  prospect  that  only  quite  modest  arrays  of  equivalent  sources  may  be  sufficient  for 
UC  and  strong  clutter  suppression  indicates  that  the  method  should  be  tested  on  data 
from  such  new  instruments  as  those  in  SERDP  project  MM-1634  and  ESTCP  MM- 
0601.  The  latter,  for  example,  involves  an  instrument  with  a  five  by  five  array  of  rather 
large  Rx  loops,  with  inherently  limited  resolution.  However,  that  may  not  be  important 
for  UC.  The  sensor  also  provides  something  of  a  blend  of  a  mono-static  data  and  the  fully 
multi-static  measurement  as  considered  in  this  study.  While  the  primary  field  is 
transmitted  in  sequence  from  different  positions  and  is  therefore  changing,  for  each 
primary  field  a  number  of  spatially  distributed  observations  are  taken.  In  the  very  least, 
we  should  investigate  ways  to  transform  data  from  both  projects  by  superposition  into  an 
equivalent  distributed  response  to  a  broad  single  excitation,  ripe  for  UC. 

5.  This  last  point  leads  to  the  most  inspiring,  but  also  most  challenging  prospect,  namely, 
generalization  of  UC  type  clutter  suppression  to  mono-static  data,  such  as  that  from 
the  EM-61,  EM-63,  GEM  instruments,  et  al.  The  formulations  would  have  much  in 
common  with  those  generated  here,  but  would  also  be  fundamentally  different  in  their 
meaning.  If  one  obtains  a  distribution  of  mono-static  signals  over  a  plane,  say,  and  then 
solves  for  an  underlying  distribution  of  equivalent  sources  over  an  arbitrarily  deeper 
plane  that  will  reproduce  these  data,  he  is  really  inferring  an  equivalent  object  below,  not 
a  boundary  condition.  That  is,  the  inferred  distribution  of  sources  doesn’t  simply  succeed 
in  reproducing  the  values  over  a  surface  for  a  single  field,  and  hence  allowing  calculation 
of  the  rest  of  the  field  (elsewhere).  Rather,  the  equivalent  source  distribution  over  (say) 
the  surface  of  the  ground  must  somehow  contain  all  the  responsiveness  of  whatever 
physical  sources  are  distributed  in  3-D  below  the  ground.  For  example,  one  might 
construct  an  array  or  continuous  distribution  over  the  ground  of  magnetic  dipoles.  Their 
responses  would  depend  on  the  primary  field  locally  striking  them,  which  might  be  made 
to  reproduce  the  responses  of  unknown  objects  below.  At  the  outset,  it  is  not  at  all  clear 
that  this  can  be  done.  Preliminary  work  done  by  our  team  using  Normalized  Surface 
Magnetic  Charge  (NSMC)  [11-13]  and  related  formulations  has  been  suggestive 
however.  Because  the  effective  suppression  of  clutter  offers  such  an  enormous  benefit. 
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this  possibility  should  be  investigated.  Then  many  well  developed  discrimination 
methods  for  mono-static  data  can  be  applied  much  more  effectively. 

6.  Not  so  well  developed  discrimination  methods  should  be  investigated  as  well.  In 

particular,  note  that  one  virtue  of  the  UC  computations  is  that  no  model  of  a  subsurface 
target  is  required.  No  laborious  optimizations  are  required  for  target  location,  orientation, 
or  properties  in  order  to  project  the  response  to  another  elevation.  One  simply  solves  by 
direct  computation  for  equivalent  source  distributions  over  an  arbitrarily  chosen  surface 
below  the  data.  This  inspires  the  notion  of  further  ‘^model-free”  processing,  leading 
beyond  signal  clean-up  to  UXO  discrimination.  We  have  investigated  such  methods 
based  on  straightforward,  robust  integral  measures  of  signal  properties,  such  as  total 
intensity  combined  with  measures  of  spatial  distribution.  Analytical  expressions, 
empirical  relations,  and  fancier  statistical  learning  machine  experiments  suggest  that  one 
can  infer  from  these  signal  measures  the  basic  target  properties  one  needs  for  dig/  no-dig 
decisions:  location  and  size  of  target.  This  kind  of  discrimination  processing,  which 
requires  no  target  model,  optimization,  or  signal  matching,  would  be  performed  on 
the  UC  signal,  which  would  be  largely  devoid  of  distortion  due  to  clutter  and  noise. 
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2.  Governing  Physics  and  Mathematical  Relations 


Magnetometry  and  electromagnetic  induction  (EMI)  sensors  operate  in  the  magneto¬ 
quasistatic  regime  (0-100’s  of  kHz).  Displacement  currents  are  negligible  [5].  Further,  in  our 
typically  shallow  UXO  discrimination  surveying,  electric  currents  in  the  soil  are  some  nine  or 
more  orders  of  magnitude  weaker  than  those  in  targets  of  interest,  so  that  the  latter  dominates  the 
fields  received  by  a  sensor.  This  is  because,  according  to  the  basic  continuity  laws  of 
electromagnetics  [5,6],  electric  fields  in  a  metallic  target  and  in  the  nearby  soil  are  similar  in 
magnitude.  However  electrical  conductivity  for  UXO  metals  is  on  the  order  of  10^  S/m,  while  that 
of  the  soil  is  likely  to  be  at  most  -  10'^  S/m  [6].  Hence  both  electric  and  displacement  currents  are 
insignificant  in  the  soil.  Both  currents  types  will  also  be  negligible  in  the  air,  which  possesses 
electrical  conductivity  that  is  much  less  than  that  of  the  weakly  conducting  soil  [6].  Altogether, 
this  means  that,  in  the  region  where  a  sensor  resides  outside  a  metallic  target.  Ampere’s  Law  in 
Maxwell’s  governing  equations  for  electromagnetics  [5,6]  becomes. 


VxH 


Vx 


UJ 


0 


(3) 


For  our  purposes  here,  the  magnetic  permeability  p  is  constant  and  can  be  factored  out  of  the 
equation.  Thus  both  the  magnetic  flux  intensity  (B  field,  in  teslas)  as  well  as  the  magnetic  field 
intensity  (H,  in  A/m)  are  irrotational.  That  in  turn  implies  that  each  can  be  expressed  by  the 
gradient  of  a  scalar  potential,  i.e. 

B  =  -V^  (4) 

Together  with  Gauss’s  law  for  magnetic  fields 

V-B  =  (5) 


this  implies  a  governing  scalar  Laplace  equation. 
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(6) 


=  -q^ 


where  physically  the  net  magnetic  charge  density  is  actually  always  zero.  While  nature  does 
not  permit  the  existence  of  isolated  magnetic  monopoles,  i.e.  positive  or  negative  points  of  q^,  in 
some  instances,  as  below,  inclusion  of  this  quantity  via  “equivalent”  sources  may  prove  to  be  a 
useful  mathematical  device.  For  example,  the  basic  physically  realizable  magnetic  source  form  is 
a  dipole.  While  such  a  structure  arises  from  electric  current  loops,  it  can  be  represented 
conveniently  in  regions  away  from  the  source  by  the  pairing  of  positive  and  negative  fictitious 
magnetic  charges.  These  produce  the  correct  fields  in  the  sense  that,  in  a  region  of  interest  away 
from  the  source,  the  fields  satisfy  the  governing  homogeneous  version  of  (6)  and  the  relevant 
boundary  conditions.  If  is  retained  as  a  source  term,  then  the  solution  of  (6)  is  [5,6] 


(7) 


where  5*0  is  surface  on  which  the  equivalent  charges  reside  and  r  is  an  observation  point. 

Overall,  our  instruments  measure  magnetic  field,  not  potential.  Given  at  least  some  vector 
components  of  B  measured  over  a  surface,  one  could  in  principle  use  the  data  as  a  boundary 
condition  on  (3)  and  (5),  with  q^  =  0,  and  then  solve  for  B  in  the  region  above.  However  this 
presupposes  rather  complete,  high  resolution  boundary  data  in  B,  particularly  if  differential  forms 
are  applied.  Finite  differencing  of  the  boundary  data  will  generally  be  severely  noise-amplifying. 
Various  integral  equation  solutions  for  B  above  a  measurement  plane  are  also  possible.  In 
general,  these  also  require  manipulation  of  vector  components  of  B  and  possibly  derivatives  or 
spatial  transforms  of  those  components  to  obtain  the  complete  information  needed. 

The  simplest  and  probably  most  robust  way  to  express  B  is  via  (4)  and  (7), 

B(r)  =  fdV  (g) 
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where  R  =  r-r',  R  =  R/|r|  and  the  equivalent  sources  do  not  reside  within  the  region  of  interest 

although  they  may  be  on  its  boundary.  This  is  an  integral  equation  formulation,  but  in  not  in 
the  sense  of  integral  equation  formulations  only  involving  fields.  A  set  of  equivalent  sources  in 
(8)  will  inherently  satisfy  the  divergence-free,  curl-free  requirements  of  our  governing  equations 
within  the  region  of  interest.  If  they  further  satisfy  a  legitimate  boundary  condition  of  that  region 
then  they  provide  all  that  one  needs  for  complete  definition  of  the  field  at  hand. 

In  this  study,  only  time  domain  (TD)  instruments  were  used,  so  that  B  above  is  actually 
B(r,^).  However,  this  is  not  significant.  One  can  express  the  B  distribution  at  any  given  time  just 
on  the  basis  of  boundary  conditions  or  source  values  at  that  same  point  in  time.  That  is,  time  does 
not  appear  explicitly  in  (3)  -  (6).  It  only  enters  implicitly  through  forcing  functions,  i.e.  boundary 
conditions  or  equivalently  through  q^.  The  time  and  space  scales  in  the  MQS  regime  are  such  that 
there  is  no  significant  delay  between  changes  through  time  on  a  boundary  or  in  a  source  and  those 
that  result  from  them,  some  distance  away.  Altogether,  noting  that  the  formulation  applies 
whether  q^  is  a  function  of  time  or  frequency,  we  assume  here  that  a  time  dependence  is  implied 
and  suppress  its  expression  in  the  equations. 

Consider  a  distribution  of  charge  q^  over  a  2-D  the  surface  So  (Figure  2),  which  in  principle 
need  not  be  flat.  As  these  are  equivalent  sources,  one  may  consider  that  they  produce  fields  as  if 
they  resided  in  infinite  homogeneous  space,  i.e.  (8)  applies  without  consideration  of  material 
properties,  discontinuities,  etc.  The  q^  could  form  a  continuous  function,  or  be  a  set  of  delta 
functions,  as  shown  and  as  implemented  computationally  in  this  study. 


I 


B  (r-) 


Figure  2.  A  source  distribution  of  equivalent  magnetic  charges  ^Bover  a  flat  surface. 
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Standard  limiting  arguments  based  on  (5)  [5,6]  produce  a  relation  stating  that  any  jump  in  the 
normal  component  of  B  across  the  surface  is  equal  to  the  magnitude  of  the  source  concentration 
on  that  surface.  Thus,  for  a  continuous  distribution 


^n(r")-5„(r')  =  (9) 

where  ^  and  ^  are  normal  field  component  values  on  opposite  sides  of  the  surface. 

The  source  concentration  at  a  point  radiates  equally  in  all  directions,  in  particular  equally 
into  the  regions  above  and  below  the  surface.  Thus,  as  suggested  in  the  figure,  BJj^)  and  BJj~) 
are  equal  in  magnitude  but  opposite  in  direction,  so  that 

=  25„(r")  (10) 


If  one  can  obtain  sufficiently  well  sampled  B^  values  over  a  surface,  this  is  very  convenient  for 
computing  the  field  elsewhere  as  (10)  is  pointwise  applicable.  That  is,  given  a  5n  value  at  any 
point  on  the  surface  immediately  allows  the  specification  of  there  without  further  solution. 
Information  on  this  single  field  component  over  So  allows  immediate  inference  of  sources  that 
are  sufficient  for  all  components  of  B  beyond  5*0.  However,  note  that  the  component  used  must 
be  the  normal  component,  and  the  sources  implied  reside  on  the  surface  immediately  adjacent  to 
the  positions  where  B^  was  obtained.  As  explained  in  the  next  section,  we  obtain  as  data  here  the 
surface  distribution  of  a  B  value  component  that  is  in  a  consistent  but  arbitrary  direction  relative 
to  the  surface  normal.  This  does  not  allow  us  to  infer  directly,  point  by  point  as  per  (10),  the 
value  of  qB  on  the  data  surface. 

To  proceed  nonetheless,  consider  an  arrangement  with  magnetic  field  data  B(r)  on  some 
measurement  surface  together  with  associated  equivalent  sources  distributed  on  an  offset 
surface  Sq,  as  shown  schematically  in  Figure  3. 
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B(r) 


Figure  3.  Charge  sheet  offset  from  observation  (data)  point. 


In  contrast  to  the  relation  in  (10)  in  which  the  charge  value  at  a  single  point  implies  the  value  of  a 


B  component  there,  here  all  charges  contribute  to  each  point  value  of  B.  If  the  charge  distribution 
is  correct  (i.e.  such  as  would  be  obtained  from  (10)  applied  on  Sq),  it  will  provide  all  components 
of  B  on  Conversely,  if  any  of  the  B  components  is  changed  over  that  must  imply  some 
change  in  the  q^.  Thus,  from  data  on  an  offset  surface  one  should  be  able  to  solve  for  the 
underlying  on  Sq  using  any  consistent  distribution  of  values  for  a  given  component  of  B. 
Readers  left  a  bit  queasy  by  this  rather  sketchy  treatment  of  uniqueness  issues  are  referred  to 
comments  elsewhere  in  this  report.  As  presented  thus  far,  these  remarks  may  at  least  serve  as  an 
informative  context  and  motivation  to  follow  the  computational  procedures  successfully  applied 
in  the  project. 

When  a  component  of  the  response  B  field  in  some  arbitrary  B^  direction  constitutes  the 
data,  (8)  provides  the  signal  as. 


(11) 


which  is  discretized  computationally  as 


(12) 
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The  source  distribution  in  (1 1)  has  been  treated  as  a  collection  of  delta  functions  in  space,  of 

magnitude  at  V-  ;  Si  is  the  signal  at  ri  and  Rji  is  the  distance  from  r.  to  ri. 

j  J  J 


Lastly,  let  us  note  an  important  distinction  regarding  the  kind  of  data  that  can  be  treated  in 
continuation  calculations.  Most  EMI  instruments  enlisted  for  UXO  discrimination  studies  or  field 
work  produce  at  least  approximately  mono-static  data.  That  is,  the  signal  is  received  at  the  same 
location  from  which  the  primary  (excitation)  field  is  transmitted.  This  means  that,  as  the  sensor  is 
moved  around  an  object  below  ground,  each  observation  represents  a  view  from  a  different 
position  of  the  target’s  response  to  excitation  from  that  position.  Thus  each  recorded  response  is 
a  sample  of  a  different  response  (“secondary”)  field.  This  is  fundamentally  different  from  the 
situation  that  we  have  considered  thus  far  and  will  consider  everywhere  below.  In  particular,  we 
assume  that  the  excitation  field  is  invariant  in  the  sense  that  it  does  not  change  as  a  result  of 
different  observation  positions  of  the  instrument.  In  traditional,  static  magnetometry  the 
excitation  field  is  the  earth’s  field.  On  the  time  and  space  scales  of  our  surveys,  this  may  be 
regarded  as  constant. 

In  our  field  studies  in  which  continuation  is  applied  (to  our  knowledge)  for  the  first  time  to 
transient  EMI  data,  the  excitation  field  is  provided  by  a  large  loop  surrounding  the  entire  site.  Its 
excitation  field  is  unchanged  as  surveyors  change  position  to  observe  the  response  field.  In  either 
this  case  or  that  of  static,  traditional  magnetometry,  around  a  given  target  all  data  consist  of 
different  samples  of  the  same  response  field,  which  may  vary  in  magnitude  and  direction  from 
point  to  point.  Thus,  a  sufficient  sampling  of  this  field  over  space  can  provide  a  legitimate 
boundary  condition  on  it,  or  would  allow  computation  a  set  of  equivalent  sources  that  would 
serve  the  same  function.  By  contrast,  a  distribution  of  mono-static  EMI  data  over  a  surface  (e.g. 
from  an  EM-61  or  GEM-3  instrument),  no  matter  how  complete  or  well  resolved,  does  not 
comprise  a  sufficient  boundary  condition  in  the  usual  sense.  One  could  not  use  such  mono-static 
data  to  compute  data  values  elsewhere,  at  least  as  described  above.  Details  of  the  fully  multi¬ 
static  observations  engaged  here  are  presented  in  the  Section  3.2. 
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2.1. 


Nature  of  the  TD  Magnetometer  Data 


The  time  domain  (TD)  GAP  Geophysics  system  [8],  developed  from  the  Sub-audio 
Magnetics  technology,  functions  by  surrounding  a  field  of  as  much  as  several  acres  with  a 
transmitting  (Tx)  loop.  Responses  from  objects  within  that  area  constitute  the  secondary  field, 
measured  via  local  roving  magnetometers.  While  the  receivers  register  a  time  dependent  signal, 
they  are  nonetheless  magnetometers  of  the  sort  used  in  static  magnetometry.  However  here, 
instead  of  an  essentially  invariant  earth’s  field  serving  as  the  excitation,  the  imposed  on-off  Tx 
loop  current  provides  a  single,  large  primary  field  that  does  not  vary  appreciably  over  the  space 
around  the  targets.  Decaying  target  responses  are  recorded  during  the  off  phase  of  the  Tx 
currents,  when  the  excitation  field  is  negligible.  Despite  its  TD  character,  the  transient  GAP  data 
has  fundamental  links  to  the  earth’s  field  in  a  manner  reminiscent  of  traditional  magnetometry 
data.  This  results  from  a  combination  of  the  receiver  type  and  the  fact  that,  while  the  Tx  current  is 
off  during  the  recording  time,  the  earth’s  field  is  not. 

The  cesium  type  magnetometer  receivers  in  the  GAP  system  respond  fundamentally  to  the 
magnitude  of  the  total  magnetic  field  that  impinges  upon  them.  That  is,  the  resonances  in  the 

receiver  respond  to  the  total  |b(^)|  =  [t)  +  {t)  +  Bl  {t)  including  the  earth’s  field. 

However,  continuation  calculations  require  at  least  partial  information  on  individual  vector 
components.  If  one  begins  with  an  equation  such  as  (8)  but  manipulates  it  so  as  to  proceed  in 
terms  of  |B|  instead  of  B  or  its  components,  the  problem  becomes  nonlinear.  In  addition  to  adding 
much  complexity,  the  optimizations  required  to  solve  for  the  would  be  vulnerable  to  non¬ 
uniqueness  problems  as  well.  In  practice,  as  explained  below,  the  inclusion  of  the  earth’s  very 
large  field  in  the  raw  data  compensates  for  added  complexity  by  allowing  some  beneficent 
approximations  [9],  which  in  turn  permit  one  to  proceed  with  continuation. 

The  perturbation  field,  b,  produced  by  the  target  of  interest  in  response  to  the  GAP  Tx  field 
is  much  smaller  than  the  earth’s  field,  Bg. 

Earth' s  field  magnitude  =  |bJ  =  B^  »  b  =  |b|  =  anomaly  magnitude  (13) 

The  signal  S  from  the  GAP  receivers  consists  essentially  of  the  magnitude  of  the  earth’s  field  plus 
the  perturbation,  minus  the  magnitude  of  the  earth’s  field  alone. 
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B  +b  -S 


(14) 


Figure  4  shows  the  vectors  involved.  This  is  only  a  schematic  diagram  in  that  the  origin  is 
substantially  farther  away  than  is  depicted,  so  that  the  two  vectors  containing  the  earth’s  field  are 
longer  and  more  parallel  to  each  other  than  is  shown.  In  any  case,  even  as  the  figure  appears,  one 
sees  geometrically  that  the  signal  consists  effectively  of  the  projection  of  the  perturbation  b  onto 
the  direction  of  the  earth’s  field. 


Figure  4.  Schematic  vector  diagram  of  the  field  impinging  on  the  receivers,  consisting  of  the 
earth’s  field  with  and  without  the  target  perturbation  b. 


Expressed  in  equations, 

|B,+b|-5,  =  7(B,+b)-(B,+b)  -  B, 

«  B^  Vl  +  2Be  -b/^e  - 

«  5,(l  +  B,.b/5^+A.oi)-  B^  «  B,.b 


(15) 


where  the  higher  order  terms  Qi.o.t.)  in  the  Taylor  series  are  easily  neglected.  With  the  data 
consisting  essentially  of  the  component  of  the  perturbation  field  in  the  earth’s  field  direction,  S  = 
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Z?e,  the  governing  equation  just  becomes  (1 1),  i.e.  numerically  (12).  This  is  the  relation  used  in  all 
the  computations  that  follow.  A  value  for  applicable  at  Aberdeen  Proving  Ground  testsite  was 

obtained  from  the  USGS  [http://geomag.usgs.gov/],  which  furnishes  detailed  data,  charts,  and 
models  easily  capable  of  producing  more  accurate  values  than  are  needed  here.  Using  the  IGRF- 
2005  model  from  that  source  for  6/15/2007,  lat  39N31  =  39.5167,  long  =  76W08  =  -76.13 
produces  components  of  the  earth  field 

=  20,149  nT 

=  -4,092  nT  (16) 

=  48,437  nT 

which,  when  translated  into  the  coordinate  system  used  in  the  analyses  (X  positive  north,  Z 
positive  up),  produces 

=  [-0.0778,0.3830,-0.9205]  (17) 


3.  Upward  Continuation 


3. 1.  Simulations 


Prior  to  field  testing,  simulations  were  performed  to  assist  in  algorithm  development;  to 
validate  the  concepts  and  computer  codes;  and  to  identify  problems  and  assess  the  usefulness  of 
the  approach,  especially  for  clutter  suppression.  A  target  and  a  nearby  clutter  item  were  assumed 
to  behave  as  anisotropic  point  magnetic  dipoles,  producing  responses  according  to  the  relation 


ItUr  \  -  I 


m,  + 


3R„R„-I 


cm  cm 


47rR: 


•in„ 


(18) 


17 


where  is  the  secondary  field  at  a  measurement  point  r„i;  the  subscript  “tm”  indicates  target-to- 
measurement  point,  while  “cm”  indicates  clutter  item-to-measurement  point;  and  nit  and  nic  are 
the  induced  target  and  clutter  dipole  moments,  respectively.  The  primary  field  that  stimulates  nit 
and  nic  does  not  appear  explicitly  in  the  expressions  because,  as  in  the  field  setting,  it  is 
essentially  uniform  and  vertical.  We  assume  that  the  clutter  item  has  an  inherently  smaller  dipole 
moment  than  the  larger  target  but  that  its  signal  obtrudes  because  it  is  shallower.  The  signal 

pattern  is  complicated  as  well  by  the  fact  that  it  is  sampled  in  a  generally  oblique  direction,  as 

per  (11).  This  can  cause  changes  in  signal  sign  from  one  side  of  the  target  to  the  other  even 
although  the  secondary  field  is  largely  vertical.  Figure  5  show  signal  contours  for  a  hypothetical 
situation  in  which  nic  ^  [0,0,10  ^1,  mt  ^  [0,0,1],  and  the  measurement  plane  is  at  Zm  =  0.2  m.  The 
target  is  beneath  the  origin  at  z  =  -0.6  m  while  the  clutter  item  is  shallower  and  is  offset,  i.e.  Yc  = 
[0.5,  0.5,  -0.2]  m.  The  contrast  between  the  top  and  bottom  plots  shows  the  asymmetries 

introduced  by  the  sampling  of  relative  to  what  one  would  see  just  by  apprehending  the 
vertical  field  component. 

Figure  6  shows  an  upward  continuation  test  for  this  situation.  The  top  plot  shows  again  the 
data  on  the  measurement  plane,  i.e.  B^  component  of  the  secondary  field  at  Zm  =  0.2  m.  The 

signal  that  would  be  obtained  at  Z  =  0.4  m  is  also  computed  directly  from  the  two  actual  sources 
(middle  plot).  The  bottom  plot  is  obtained  by  first  solving  for  a  plane  of  equivalent  sources  at 
Zq  =  0  using  (12),  based  on  the  data  S  at  Zm  shown  in  the  top  plot.  This  solution  is  then  used 
to  compute  the  field  at  Z  =  0.4  m.  Clearly  the  upward  continuation  solution  does  well  in 
predicting  the  signal  at  a  higher  elevation.  Equally  important,  both  actual  and  continuation  plots 
show  a  notable  diminution  of  the  prominence  of  the  signal  from  the  shallow  clutter  item. 
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Figure  5.  Simulated  responses  of  a  target  and  clutter  item  combination,  when  is  just  z 
(top),  and  when  is  [-0.0778,  0.3830,  -0.9205]  as  per  the  APG  site  (bottom). 
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Figure  6.  Same  case  as  in  Figure  5,  showing  component  of  ff .  Top:  Data  on 

measurement  plane  (Zm  =  0.2  m).  Middle:  actual  signal  that  would  be  obtained  at 
a  higher  elevation  (Z  =  0.4  m).  Bottom:  Continuation  signal  computed  for  Z  =  0.4 
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In  these  computations  the  dimensions  of  So  were  the  same  as  those  of  (see  plots). 

Spatial  increments  between  points  in  the  square  arrays  of  data  and  of  sources  were  the  same,  0.15 
m,  being  only  slightly  less  than  the  separation  between  So  and  While  this  discretization  is 
rather  coarse,  it  still  works  well  and  is  chosen  purposely  to  avoid  ill-conditioning.  With  a 
condition  number  of  ~  75  one  would  expect  distortions  in  when  noise  at  the  level  shown  in 
Figure  7,  top,  is  added  (See  Section  5).  Here  RelNoiseMag  =  0.25,  where  that  parameter  implies 
added  white  noise  between  +! -RelNoiseMag! 2  times  the  signal  maximum.  While  distortions  in 
may  indeed  occur,  the  bottom  plot  shows  that  to  a  remarkable  extent  their  effects  are  integrated 
out  in  the  course  of  the  upward  continuation:  resemblance  between  magnitude  and  shape  of  the 
signal  distributions  in  the  middle  and  bottom  plots  is  clear.  Thus  the  q^  solution  and  continuation 
process  together  appear  quite  robust  in  the  fact  of  substantial  noise  in  the  data.  Both  random  noise 
as  well  as  physical  clutter  signals  are  suppressed  by  the  process.  Figure  8  pursues  this,  with 
somewhat  more  moderate  noise  {RelNoiseMag  =  0.1)  and  nit  =  [0,1,1]  for  a  more  complicated 
response  pattern.  The  system  still  works  equally  well. 

Of  course,  in  all  this  we  are  primarily  interested  in  isolating  the  signal  from  the  target  itself 
Figure  9  illustrates  the  difference  between  a  signal  on  the  data  plane  from  the  target  alone, 
compared  to  that  from  the  combination  of  target,  shallow  clutter  item,  and  random  noise.  When 
upwards  continuation  processing  is  applied  to  the  noisy  and  cluttered  data,  the  results  in  Figure  10 
are  obtained.  Continuation  has  done  a  good  job  of  predicting  the  signal  from  target  and  clutter 
item  at  considerable  elevation.  While  the  magnitude  from  this  combination  is  still  slightly  high 
relative  to  that  from  the  target  alone,  the  differences  are  minimal.  For  all  practical  purposes,  the 
method  shows  promise  for  isolating  target  contributions  in  a  cluttered,  noisy  signal. 

As  mentioned  above,  noise  in  the  data  can  distort  the  q^  solution,  but  in  the  examples  here 
that  did  not  ultimately  impact  the  quality  of  upward  continuation  significantly.  The  distortions  are 
still  worrisome,  however,  in  that  under  some  other  circumstances  they  could  carry  through  to  the 
continued  signal.  Further,  if  one  wishes  to  examine  the  q^  solution  itself,  the  distortions  could  be 
problematical.  The  distortions  in  q^  can  be  especially  prominent  when  one  chooses  geometrical 
parameters  different  from  those  above,  most  notably  deeper  So,  as  for  the  downward 
computations  that  are  one  investigatory  focus  of  this  project.  Figure  1 1  shows  the  q^  solution 
obtained  in  the  problem  we  have  been  examining,  with  and  without  noise  added  to  the  data 
{RelNoiseMag  =  0.1).  In  that  the  continuation  results  are  good,  the  essential  source  distribution 
is  clearly  buried  within  the  heart  of  the  distorted  solution,  whose  rapid  oscillations  are  ultimately 
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integrated  out  in  the  continuation  integration  (12).  Nevertheless,  the  lower  plot  is  not  a  pretty 
picture.  Especially  if  distorted  further,  its  pattern  would  be  difficult  to  interpret  and  its  errors 
could  threaten  the  quality  of  upward  inferences.  Section  5  discusses  the  nature,  origins  and 
control  of  such  distortions. 
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Figure  7.  Same  case  as  above  but  with  RelNoiseMag  =  0.25.  Top:  Data  on  at  Zm  =  0.2 
m.  Middle:  Actual  signal  from  targets,  without  noise,  at  Z  =  0.6  m.  Bottom: 
Signal  at  Z  =  0.6  m  obtained  by  continuation  from  based  on  data  in  top  plot. 
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Figure  8.  Same  as  previous  plot,  but  with  RelNoise  =  0.1  and  mt  =  [0,1,1]. 
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Figure  9.  Same  in  previous  figure  but  showing  difference  between  data  with  noise  and  also 
clutter  item  (top)  and  from  target  only  (bottom). 
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Actual  Be  comp  on  continuation  plane  Z  =  1.2,  no  clutter/noise 
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Actual  Be  comp  on  continuation  plane  Z  =  1.2,  w/  clutter,  no  noise 
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Figure  10.  Same  parameters  as  in  Figure  9,  showing  little  difference  at  elevation  between 

theoretical  target-only  signal  (top),  that  from  the  target  plus  clutter  (middle),  and 
that  from  continuation  using  target  +  clutter  +  noise  on  the  data  plane  (bottom). 
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Figure  1 1 .  The  solution  for  the  same  case  but  with  RelNoiseMag  =  0  (top)  and 
RelNoiseMag  =  O.l  (bottom). 
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3.2. 


Field  Tests 


During  June  2007  the  team  made  measurements  with  the  GAP  Geophysics  survey  system 
[8]  at  the  UXO  Standardized  Test  Site  at  Aberdeen  Proving  Ground  (APG).  A  section  of  field 
approximately  40  m  square  was  surrounded  by  the  transmitter  coil  and  measurements  were  made 
within  it  along  lines  over  an  area  approximately  4  m  by  9  m.  The  wheeled  rig  proceeded  along 
lines  guided  visually  by  taught  strings.  Figure  12  shows  the  rig  proceeding  along  a  line  and  the 
pattern  of  receiver  positions  along  all  the  lines.  The  rig  transported  four  receivers  spaced  38  cm 
apart.  Lines  were  spaced  so  that  the  sensor  tracks  were  interleaved,  producing  on  average  an 
approximately  19  cm  data  spacing  in  the  across  track  direction.  Data  along  track  were  taken 
approximately  every  0.3  sec  which,  with  the  slow  rate  of  progress  maintained,  meant  a  1  to  2  cm 
data  density  in  the  along  track  direction. 


Figure  12.  Left:  survey  rig  proceeding  along  a  survey  line  across  the  edge  of  the  sand  pit  in 
which  targets  were  buried.  Right:  the  paths  of  the  individual  receivers. 


28 


Targets  were  buried  in  about  the  center  of  the  sand  pit  (hole  in  Figure  12,  close  up  in  Figure 
13).  Primarily,  a  105  mm  projectile  was  used  in  various  horizontal,  vertical,  and  inclined 
orientations. 


Figure  13.  Hole  dug  in  sand  pit  with  105  mm  projectile  placed  horizontally  at  bottom  and 
clutter  item  (shotput)  on  surface,  with  meter  measure. 


Positioning  data  was  provided  by  a  GPS  antenna  attached  to  the  rig,  indicated  in  Figure  12 
and  Figure  14.  This  was  elevated  above  the  receiver  platform  to  avoid  electronic  interference.  It 
was  lowered  from  its  initially  high  position  (Figure  12)  to  that  in  Figure  14  in  order  to  minimize 
positional  error  due  to  tilt  of  the  platform  on  the  uneven  ground.  Examination  of  the  data  along 
the  lines  suggests  a  residual  random  error  from  platform  motion  of  about  1  or  2  cm. 

As  the  data  density  along  and  across  track  are  quite  different,  a  specialized  interpolation 
algorithm  was  devised  to  produce  2-D  grids  of  data  with  uniform  spacing  in  both  horizontal 
directions.  While  data  produced  by  the  algorithm  should  not  possess  any  greater  positional  error 
than  the  original,  it  will  not  be  improved  either.  This  should  be  bom  in  mind  when  interpreting 
plots  such  as  the  example  signal  in  Figure  15.  Drawing  contours  (i.e.  via  Matlab)  can  produce  an 
illusion  of  isotropic  resolution,  but  examination  of  the  horizontal  vs  vertical  patterns  in  the  plot 
reveals  the  remains  of  the  originally  unequal  data  densities. 
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Figure  14.  Rig  passing  over  target  region,  with  target  hole  covered  by  plywood  and  one  (top) 
or  three  (bottom)  shotputs  placed  nearby  to  serve  as  shallow  clutter  items. 
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Figure  15.  Data  from  lines  such  as  those  in  Figure  12  interpolated  onto  a  square  grid  of 
points  at  10  cm  x-y  spacing. 


The  spatial  distributions  of  data  such  as  that  in  Figure  15  pertain  to  a  particular  point  in 
time.  At  any  given  point  in  space,  the  sensor’s  receivers  see  a  time  pattern  of  magnetic  fields 
illustrated  in  a  plot  of  field  data  in  Figure  16,  top,  from  a  relevant  portion  of  the  on-off  Tx  cycle. 
The  earth’s  field  contribution  has  been  subtracted.  When  the  voltage  in  the  Tx  cable  is  switched 
on,  before  time  zero,  the  current  and  consequent  Tx  field  rise  to  an  essentially  flat,  high  value, 
which  gradually  saturates  the  target.  When  the  Tx  voltage  is  shut  off  at  time  zero,  the  primary 
field  decays  rapidly  until  only  the  secondary  field  from  the  target  remains  at  the  receiver.  The 
log-log  decay  of  the  latter  is  shown  in  Figure  16,  bottom.  Note  that  the  magnitude  of  the  earth’s 
field  is  ~  50,000  nT.  The  magnitudes  of  signals  shown  in  Figure  16  validate  the  assumption  that 
the  signals  of  interest  -  the  perturbation  quantity  b  above  -  are  indeed  much,  much  smaller  than 
the  ambient  earth’s  field.  The  first  time  gate  in  the  plot  and  in  data  analyzed  is  chosen  to  be  late 
enough  to  avoid  distortions  associated  with  the  lingering  decline  of  the  primary  field  in  earliest 
time  after  zero.  In  all  analyses  here,  which  pertain  to  spatial  distributions  of  signals  at  a  given 
time,  it  does  not  matter  much  in  principle  what  time  point  is  chosen  during  the  decay  in  Figure 
16,  bottom.  In  an  actual  survey  one  might  wish  to  examine  later  time  points  at  which  responses  of 
smaller  clutter  items  might  well  have  diminished,  relative  to  that  from  larger  items  of  interest.  To 
retain  focus  on  clutter  vs  target  signals,  it  behooves  us  here  to  use  early  time  points  where  both 
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signals  are  at  their  clearest.  Therefore  in  all  treatments  of  field  data  below  we  choose  a  time  gate 
at  0.625  ms  after  time  zero. 


Full  Waveform 


Time  (s) 


Figure  16.  Data  taken  over  the  peak  of  the  target  signal  for  a  105  mm  projectile  buried  73  cm 
below  the  sensor.  Top:  In  linear  time,  full  waveform  including  Tx  field  prior  to 
time  zero,  rapid  decline  after  shutoff,  followed  by  slow  decay  associated  with 
target.  Bottom:  Decay  of  the  target,  isolated,  relative  to  log(t). 


The  level  of  inherent  noise  in  the  data  might  be  best  apprehended  from  the  plots  of 
individual  receiver  outputs  along  the  survey  lines  (Figure  17)  rather  than  by  examination  of 
interpolated  values  in  Figure  15.  Not  surprisingly,  the  relative  magnitude  of  the  noise  is  more  or 
less  inversely  proportional  to  the  overall  signal  strength.  In  general,  relative  to  peak  data  values 
(Figure  17,  top),  the  SNR  is  good. 
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X  =  2.37 


X  =  1  .43 


Figure  17.  Data  along  survey  lines  over  a  vertical  buried  105  mm  UXO  for  the  same  case 
shown  in  Figure  15,  illustrating  the  level  of  noise. 
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For  a  first  look  at  upward  continuation  of  field  data,  consider  the  case  of  a  vertical  105  mm 
projectile  buried  initially  73  cm  to  its  shallowest  point  below  the  sensor  (Figure  18).  In  all  cases, 
because  it  was  impractical  to  redo  surveys  with  the  sensor  raised  to  a  higher  elevation,  instead  the 
same  thing  was  accomplished  by  placing  all  targets  deeper  (soil  response  was  not  significant).  At 
the  initial  distance  between  sensor  and  target,  some  of  the  apparent  lateral  (x-direction)  spread  of 
the  signal  in  the  top  plot  is  caused  by  the  rather  large  separation  of  the  survey  lines  in  that 
direction  and  interpolation  of  values  in  between.  In  any  case,  when  the  sensor-target  distance  is 
increased  to  96  cm,  the  response  is  more  clearly  centered  and  also  more  broadly  spread  (middle 
plot).  Continuation  computations  for  the  96  cm  separation,  based  on  the  initial  73  cm  data,  are 
shown  in  the  bottom  plot.  Continuation  catches  all  the  essentials  of  the  actual  elevated 
measurements  in  both  magnitude  and  spread  of  the  signal,  in  the  process  filtering  out  a  good  deal 
of  random  irregularity  in  the  signal. 

The  quality  and  resolution  of  the  data  in  this  example  and  our  ability  to  institute  only 
modest  changes  in  sensor  elevation  above  the  target  limit  the  appeal  of  this  example  as  a 
demonstration  of  the  technique  and  its  benefits.  Computations  further  below  exploit  the  real 
appeal  of  the  method,  namely,  its  ability  to  project  data  to  separations  that  could  not  be  attained  in 
the  physical  tests  with  consequent  substantial  benefits  of  clarity.  Remaining  for  the  moment  with 
more  modest  changes  in  separation  for  the  sake  of  validation,  some  results  involving  clutter  are 
illuminating.  Figure  19  shows  a  case  in  which  a  single  shotput  was  placed  at  ground  level  to  serve 
as  a  clutter  item.  At  the  initial  elevation,  both  clutter  and  UXO  are  approximately  equal  in 
magnitude.  When  the  sensor-target  distance  is  increased  (middle  plot)  signals  from  both  sources 
decline,  but  the  strength  of  the  clutter  signal  fades  relative  to  that  from  the  UXO.  This  is 
reproduced  in  the  signal  obtained  by  upward  continuation  (bottom  plot).  The  projected 
magnitudes  are  approximately  correct,  while  the  patterns  are  actually  cleaned  up.  Continuation 
calculations  will  reproduce  most  of  the  signal  irregularities  in  the  original  data  (see  next 
example).  However,  upward  continuation  will  tend  to  iron  these  wrinkles  out,  for  good 
physical/mathematical  reasons.  At  the  same  time,  actually  obtaining  signals  with  the  sensor  at 
greater  distances  overwhelms  this  physical  effect,  as  target  signals  become  weaker  relative  to  any 
irreducible  background  noise/  clutter. 
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Figure  18.  Signal  from  a  vertical  105  mm  UXO.  Top:  Data  at  =  73  cm.  Middle:  Data  at  <7  = 
96  cm.  Bottom:  projection  from  data  at  73  cm  to  signal  at  96  cm. 
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Figure  19.  Single  shotput  clutter  case  with  105  mm  projectile  buried  vertical  nose  down 

initially  at  73  cm  from  sensor  to  tail.  Top:  Data  at  initial  sensor  distance.  Middle: 
Data  at  23  cm  higher  sensor  distance  from  both  targets.  Bottom:  Computational 
continuation  23  upwards  of  data  in  top  plot. 
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An  example  with  two  clutter  items  (Figure  20)  shows  even  more  pronounced  reduction  of 
clutter  signals  with  increased  distance  from  the  sensor.  The  clutter  arrangement  includes  the  right¬ 
most  two  shotputs  in  Figure  14,  bottom.  Despite  the  fact  that  the  identical  clutter  items  were  both 
the  same  distance  below  the  sensor  platform,  the  magnitude  of  their  signals  is  somewhat  different 
because  of  their  differing  relations  to  the  survey  lines  and  sensor  paths.  The  shotput  near  the 
bottom  of  each  plot  was  located  somewhere  between  the  sensor  lines  such  that  interpolation 
between  the  lines  could  not  catch  its  true  peak.  In  any  case,  in  the  upper  plots  the  greatest  clutter 
signal  magnitude  is  on  the  order  of  50%  greater  than  that  of  the  UXO.  At  the  greater  sensor 
distance  (bottom  plots)  the  situation  is  reversed,  with  the  UXO  signal  now  roughly  twice  that  of 
the  clutter.  Note  again  that,  in  addition  to  reproducing  this  effect,  upward  continuation  has  also 
cleaned  up  the  random  noise  that  actually  increases  in  the  physical  data  taken  from  greater 
distance. 

Pursuing  this,  we  use  the  algorithm  to  continue  the  data  to  greater  sensor  distances  than 
were  possible  to  produce  in  the  field.  This  exploits  the  greatest  benefits  of  upward  continuation. 
Figure  21  shows  continuation  of  the  original  data  to  two  higher  levels,  in  the  course  of  which  the 
two  clutter  signals  fade  increasingly.  Ultimately  it  is  clear  that  most  of  the  signal  is  due  to  the 
target  of  interest,  the  UXO  (Figure  21,  bottom). 
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Figure  20.  Case  with  clutter  produced  by  two  shotputs  and  105  mm  projectile  at  center.  Left: 
data.  Right:  Continuation  results.  Top:  Initial  elevation.  Bottom:  23  cm  greater 
distance  to  sensor  for  all  targets. 
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Continuation  Be  comp  at  Z  =  0.5S,  w/  all  cluttor/noiso 


Continuation  Be  comp  at  Z  =  0.81,  w/  all  clutter/noise 
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Figure  21.  Computed  signals  for  105  mm  UXO  plus  two  surface  clutter  items,  with  data 
continued  beyond  the  35  cm  elevation  of  the  previous  figure,  to  58  cm,  81  cm, 
and  127  cm  (top,  middle,  bottom). 
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Before  proceeding  to  the  next  field  test,  we  return  briefly  to  a  modeling  exercise  to  assess 
better  what  one  might  expect  in  more  problematical  cases.  The  cases  above  treat  only  instances 
in  which  distinct  clutter  appears  at  positions  offset  from  the  UXO  response.  While  this  was  useful 
for  ascertaining  the  relative  strengths  of  the  clutter  and  target  signals,  it  does  not  treat  the 
situation  of  greater  interest,  namely,  when  the  clutter  is  superposed  intrusively  on  the  UXO 
signal.  Figure  22  shows  one  rendering  of  field  data  for  a  case  in  which  the  clutter  was  placed 
directly  above  the  center  point  of  an  inclined  UXO.  The  combination  of  target  inclination  and  the 

applicable  Bg  at  APG  causes  the  clutter  to  appear  off  center  on  the  target  signal  distribution.  In 

any  case,  the  clutter  produces  an  impressive  spike  on  the  smoother  UXO  response.  Can  the 
method  here  really  help  in  this  kind  of  situation? 


Figure  22.  Field  data  from  a  case  in  which  a  shotput  is  placed  on  the  surface  above  a  buried 
105  mm  projectile  inclined  at  45®. 


For  visual  simplicity  we  assume  in  a  simulation  that  =  -  z  and  the  target  is  a  vertical 

dipole  which  by  itself  will  produce  a  radially  symmetric  signal.  A  horizontally  offset  clutter  item 
with  5%  of  the  target’s  dipole  moment  is  placed  at  6  cm  depth,  vs  60  cm  for  the  UXO.  This 
produces  a  bump  on  the  target  signal  on  the  measurement  plane  at  20  cm  elevation  (Figure  23, 
top)  comparable  to  that  in  field  data  in  Figure  22.  At  a  height  of  1.2  m  the  continuation  signal 
catches  the  main  character  of  the  signal  sought  from  the  target  alone.  It  has  tracked  its  decline  in 
magnitude  from  -  0.3  to  -  0.03,  with  similar  spread  of  the  overall  pattern.  These  results  suggest 
that  even  strong  clutter  directly  within  the  main  target  response  can  be  reduced  by  our  method. 
Bolstered  by  this  simulation  result,  we  proceed  to  a  comparable  case  in  the  field. 
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Figure  23.  Simulation  case  in  which  a  clutter  item  is  located  such  that  its  signal  intrudes 

upon  that  of  the  UXO.  Top:  Data  plane  Sq  at  Zm  =  0.2  m  with  prominent  clutter. 
Middle:  Continuation  signal  from  UXO  plus  clutter  at  Z  =  1.2  m.  Bottom:  Target 
only  signal,  also  at  Z  =  1.2  m. 
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Figure  24  shows  results  for  a  case  involving  a  105  mm  projectile  buried  at  a  45  deg 
inclination,  with  and  without  a  shotput  placed  directly  above.  The  addition  of  the  shotput  (upper 
right)  to  the  scenario  with  only  the  UXO  (upper  left)  produces  a  spike  in  the  signal  about  two  and 
a  half  times  the  signal  magnitude  from  the  UXO  alone.  The  clutter  item  changes  the  signal 
distribution  very  markedly.  However  upward  continuation  of  the  data  including  the  obtrusive 
clutter  catches  the  main  character  of  data  we  should  expect  at  elevation  from  the  UXO  alone. 
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Figure  24.  Field  case  with  105  mm  at  45  deg  inclination  and  =  -81.5  cm,  with  and  without 
shotput  directly  above.  Top:  Data  on  5^0  at  15  cm  elevation.  Bottom:  Continuation 
calculations  for  1  m  higher  elevation.  Left:  UXO  alone  (no  shotput).  Right:  UXO 
with  shotput  above  it. 
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4.  The  Downward  Computation  Issue 


Downward  continuation  can  mean  different  things  to  different  people;  however,  all 
employing  it  seek  in  some  way  to  extend  above-surface  information  to  infer  the  subsurface 
situation.  At  least  from  the  point  of  view  of  consistency  in  terminology,  the  most  logical  meaning 
would  be  use  of  [above-]  surface  information  to  calculate  subsurface  fields  of  the  same  kind  that 
were  measured.  Additionally,  we  take  the  term  to  mean  here  true  continuation  involving  fully 
multi-static  data,  as  opposed  to  projection  to  some  other  level  of  mono-static  data:  A  single, 
spatially  widespread  response  field  is  continued.  In  this  view,  detailed  surface  information  is, 
after  all,  a  kind  of  boundary  condition.  One  can  consider  the  boundary  to  delimit  either  an  upper 
or  lower  half  space.  While  here  the  field  values  are  not  calculated  directly  using  integral  or 
differential  equations,  the  computational  framework  is  arguably  equivalent  via  the  intermediate 
device  of  equivalent  sources.  Information  at  one  level  is  used  to  solve  for  distributions  of 
equivalent  sources  over  deeper  surfaces.  This  implies  directly  a  distribution  of  the  field 
component  normal  to  the  deeper  surface,  and  thence  all  field  components.  In  any  case,  by 
computing  source  distributions  over  deeper  and  deeper  surfaces  one  might  seek  to  find  a  depth 
where  the  source  concentration  is  greatest,  taking  this  to  be  the  location  of  the  object  of  interest. 
The  motivation  can  be  just  to  locate  the  target,  or  to  infer  some  of  its  properties  by  examining  the 
focal  source  distribution. 

Aside  from  details  of  computational  strategy  and  their  various  pro’s  and  con’s,  all 
downward  continuation  strategies  are  challenged  by  the  same  issues:  1)  Separation  of  the  data 
plane  and  the  deeper  source  (continuation)  plane  tends  to  produce  ill-conditioned  systems;  and  2) 
proceeding  deeper  with  the  calculations  is  inherently  a  clutter/noise  amplifying  process.  Under 
the  conditions  with  which  one  must  contend  here  using  total  field  magnetometers,  the  source 
plane  must  be  displaced  to  some  position  below  the  data  plane  even  in  upward  continuation.  Thus 
ill-conditioning  is  a  potential  issue  there  as  well,  even  in  the  absence  of  any  more  ambitious 
lowering  of  the  source  plane.  The  second  issue  is  limited  to  downward  continuation.  The  results 
above  show  smoothing  of  surface  noise  as  the  signal  is  continued  upward.  It  is  not  surprising  then 
that,  proceeding  in  the  opposite  direction,  one  should  confront  the  opposite  phenomenon. 
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To  make  the  problem  palpable,  the  following  example  considers  noiseless  data  with  a 
single  source  of  signal  at  a  specified  depth  (Z  =  -20  cm).  There  is  no  loss  of  generality  in  that 
noise  caused  by  perturbations  of  some  kind  at  the  surface  would  simply  serve  as  a  broader 
distribution  particular  sources,  at  a  shallower  depth.  Studying  the  effects  of  a  single  source  will 
suffice  to  make  the  point.  Also,  deleting  noise  beyond  what  is  caused  by  the  discreet  object 
avoids  confounding  the  noise/clutter  amplification  issue  with  ill-conditioning’s  similar  effects. 

For  simplicity  consider  the  earth’s  field  to  be  vertical,  and  the  dipole  moment  of  the  object  at  20 
cm  depth  to  be  m  =  [0  0  1].  Using  data  from  a  plane  that  is  20  cm  above  the  “ground  surface”( 
So  at  Zq  =  0),  with  data  and  source  points  on  a  grid  with  10  cm  spacing,  one  obtains  the  source 
distribution  in  Figure  25,  top.  Lowering  Sq  to  -20  cm,  i.e.  to  the  same  level  as  the  concentrated 
dipole  source,  produces  the  middle  plot.  This  is  as  concentrated  a  source  distribution  as  the 
particular  computational  parameters  will  permit. 

The  point  of  this  exercise  is  apparent  in  the  bottom  plot  of  Figure  25.  Here  5*0  has  been 
lowered  another  20  cm  to  Zq  =  -40  cm.  The  computed  values  oscillate  wildly  and  spread  out 
spatially.  Both  the  limits  of  the  oscillations  and  any  smoothing  of  the  distribution  comprise  values 
that  are  orders  of  magnitude  higher  than  the  most  concentrated  sources  at  the  depth  of  the  object 
(middle  plot).  The  system  is  demonstrating  that  much  more  in  the  way  of  sources  is  required  at 
more  distant  regions  to  produce  the  effect  of  the  nearer  object.  Now,  if  the  only  purpose  were  to 
isolate  the  object,  one  could  perform  these  computations  but  then  simply  discard  the  results  in 
Figure  25,  bottom,  satisfied  that  the  middle  plot  represented  something  of  a  “focal”  point. 
However  the  problem  is  clear  if  one  considers  that  the  hypothetical  object  at  20  cm  depth  might 
be  simply  a  piece  of  clutter  and  the  object  of  real  interest  deeper.  Clearly  the  kind  of  source 
distribution  in  Figure  25,  bottom,  would  obscure  any  more  reasonable  distribution  associated  with 
a  discrete  target  at  comparable  or  greater  depth.  In  the  course  of  downward  continuation,  the 
computational  system  will  in  fact  tend  to  spread  and  amplify  greatly  any  near-surface 
perturbations  or  effective  sources  in  the  data,  independent  of  (similar  looking)  ill-conditioning 
effects. 

These  results  argue  strongly  against  pursuing  downward  continuation.  At  the  same  time, 
one  can  apply  special  measures  that  may  help,  and  which  in  any  case  are  likely  to  be  required  by 
ill-conditioning,  to  which  we  pass  in  the  next  section. 
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Figure  25.  Solutions  for  obtained  in  a  hypothetical  example  when  the  target  is  at  20  cm 
depth.  The  location  of  the  source  plane  Sq  is  at  depths  Zq  =  0  (top),  Zq  =  -20  cm 
(middle),  and  Zq  =  -40  cm  (bottom). 
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5.  Computational  Control  of  Solutions  for 
Distributions  of  Sources  at  Depth 


5. 1.  Uniqueness  and  Invertability  for  Sources 


For  both  the  upward  and  downward  continuation  strategies  pursued  here,  one  needs  to 
compute  planar  distributions  of  sources  at  some  level  below  that  where  data  were  obtained. 
Especially  if  one  wishes  to  gain  some  view  of  a  source  concentration  around  the  actual  target 
location,  it  is  tempting  just  to  measure  the  secondary  field  over  a  surface  at  a  chosen  elevation 
and  then,  using  the  equivalent  of  (12),  to  solve  for  a  distribution  of  sources  over  an  entire  volume 
below.  Ideally  the  distribution  of  equivalent  sources  would  concentrate  wherever  a  significant 
physical  source  was  located.  However  this  problem  is  inherently  ill-posed.  To  see  this,  consider  a 
simple  case  in  which  the  source  is  a  unit  magnitude  point  charge  at  the  origin.  On  the  basis  of  (8) 
only  a  radial  field  component  is  produced,  namely  Bj.  =  HAn/.  Using  the  formulation  above  but 
applied  to  curved  instead  of  fiat  surface,  one  can  produce  this  field  outside  a  spherical  surface 
with  some  radius  ri  by  assuming  a  distribution  of  equivalent  magnetic  charge  over  the  surface  of 
the  sphere. 


By  the  same  reasoning,  however,  one  can  also  produce  the  same  field  from  a  distribution  of 
charges  over  a  different  sphere  about  the  origin,  e.g.  with  radius  r2>  n. 


1 

271^2^ 


5^(r>r2) 


1 

471 


(20) 


Thus  at  any  r>r2>  r\,  both  distributions  of  sources  produce  exactly  the  same  field. 
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If  the  sources  are  distributed  over  a  surface  instead  of  a  volume,  however,  a  unique  set  of 
sources  is  implied  by  the  fields  observed  over  that  surface  (see  (10)  and  discussion  surrounding 
it).  One  may  specify  sources,  point  by  point,  just  based  on  direct  correspondence  to  field 
boundary  conditions  over  the  surface.  Our  computations  are  a  bit  more  indirect  than  point  by 
point  correspondence  between  2Hz  and  in  that  the  observed  field  values  on  Sm  are  offset 
spatially  from  the  surface  So  on  which  the  sources  lies.  However  the  same  logic  prevails:  a  given 
field  boundary  condition  on  So  will  produce  specific  field  values  in  the  region  above  it.  Any 
alteration  of  those  boundary  values  will  inevitably  produce  a  change  in  field  values  somewhere  in 
the  region  above,  i.e.  where  Sm  is.  Therefore  only  a  single  boundary  condition  on  So  and  single 
corresponding  source  distribution  on  it  will  produce  interior  field  values  that  match  those 
observed  over  the  entirety  of  Sm. 

We  take  the  gist  of  the  discussion  above  to  imply  that  one  may  assume  that  a  unique  set  of 
sources  is  implied  over  a  plane  by  our  observations  of  some  component  of  the  field  over  Sm-  In 
fact,  some  issues  of  invertability  and  uniqueness,  necessity  and  sufficiency  of  conditions  are  not 
addressed  in  complete  depth  in  the  foregoing  discussion.  While  reserving  full  pursuit  of  those 
issues  for  follow-on  work,  we  take  the  reasoning  above  together  with  computational  experience  to 
provide  confidence  that  the  numerical  system  (1 1)-(12)  is  inherently  invertible. 


5.2.  Ill-conditioning  and  Its  Control 


Unfortunately,  while  this  is  reassuring,  we  have  not  entirely  escaped  the  issue  of  ill- 
posedness  in  the  form  of  its  numerical  first  cousin,  ill-conditioning.  The  problem  is  that,  as  the 
surface  Sq  is  located  at  increasingly  greater  depths  relative  to  the  observed  data,  the  system 
becomes  increasingly  ill-conditioned.  To  understand  this,  consider  the  following  relation 


2  0 

J  ax  I 

0 


(21) 


This  might  apply  to  the  EMI  response  of  some  very  elongated  object,  when  the  coordinate  system 
is  aligned  with  the  axial  and  transverse  axes  of  the  object.  The  axial  and  transverse  sensor  signals 
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on  the  right  hand  side  are  related  by  the  matrix  diagonal  values  to  the  impinging  primary  fields 
along  the  those  same  axes.  Suppose  we  observed  the  signals  on  the  right  hand  side  and  wished  to 
infer  the  primary  fields  hitting  the  target  that  had  produced  them.  In  this  instance,  with  real  data, 
the  problem  would  be  quite  difficult  even  though  the  matrix  is  immediately  invertible  by 
inspection.  This  is  because,  in  our  example,  the  transverse  response  is  inherently  so  weak.  Given 

comparable  values  of  ,  the  transverse  signal  will  be  much  weaker  than  . 

may  likely  be  predominately  noise.  This  will  produce  a  greatly  distorted  value  of  when 
is  divided  by  lO  "^. 

Basically  the  matrix  in  (21),  while  perfectly  legitimate  and  consistent,  is  telling  us  not  to  try 
to  solve  for  the  second  unknown  because  its  effects  are  much  too  weakly  represented  in  the 
system  to  be  clearly  distinguishable.  The  difficulty  is  also  evident  if  one  simply  rotates  the 
coordinate  system  so  that  the  matrix  becomes 

r  1.0  0.99991 

A  =  (22) 

0.9999  1.0 


Here  the  pathological  nature  of  the  system  is  clear  in  another  way:  the  original  system  is 
equivalent  to  having  two  nearly  identical  equations,  a  bad  prospect  for  solving  for  two  unknowns. 
For  all  practical  purposes,  the  system  really  only  contains  one  independent  statement. 

The  degree  of  ill-conditioning  in  a  system  is  typically  quantified  in  terms  of  condition 
number  C,  which  equals  the  magnitude  of  the  ratio  of  the  largest  to  the  smallest  eigenvalues  of 
the  matrix.  In  this  example  C  =  2x  lO"^,  as  the  matrix  in  (21)  in  fact  contains  the  eigenvalues  on 
the  diagonal.  A  rule  of  thumb  is  that  a  condition  number  on  the  order  of  10^  will  amplify 
unacceptably  noise  that  is  on  the  order  of  10’^  or  larger  relative  to  significant  values  of  the  data. 
That  is,  for  stable,  accurate  solutions,  one  requires  an  SNR  at  least  on  the  order  of  C,  ideally  an 
order  of  magnitude  better.  This  is  a  challenging  requirement.  For  problems  such  as  those 
considered  here,  C  values  greater  than  100  appear  very  easily,  whereas  an  SNR  greater  than  100 
is  normally  very  hard  to  come  by. 

Referring  back  to  Figure  3,  one  can  see  how  ill-conditioning  may  arise  in  applying  (12) 
here.  If  two  source  points  on  So  are  too  close  together,  they  will  produce  very  nearly  the  same 
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equation  in  the  governing  matrix:  each  has  nearly  the  same  influence  on  the  data.  For  a  given 
separation  of  source  points,  moving  So  further  away  from  Sm  has  the  same  effect  as  moving  the 
source  point  closer  together  in  terms  of  our  ability  to  distinguish  their  influence.  By  analogy,  if 
each  source  in  a  group  of  such  source  points  emits  a  glow  of  light,  when  they  are  moved  far 
enough  away  from  the  observing  surface  they  will  simply  appear  as  a  single  blob  of  glow  with  no 
distinguishable  effects  attributable  to  the  individual  sources.  Similarly,  placing  the  observation 
points  on  S^  too  close  together  will  have  the  same  effect.  Overall,  trying  to  increase  resolution  by 
using  finer  divisions  on  either  Sq  or  will  exacerbate  ill-conditioning  in  a  highly  nonlinear 
fashion. 

In  the  field  cases  considered  in  this  report,  care  was  taken  in  terms  of  separation  of  So,  Sm, 
source  points,  and  observation  points  so  that  ill-conditioning  was  not  a  problem.  However,  only 
small  changes  in  these  spacings  would  be  required  to  produce  substantially  worse  conditioning. 
Overall,  in  general  application  of  the  methods  of  upward  continuation  proposed  here  over  any 
variety  of  situations,  the  appearance  of  ill-conditioning  problems  is  virtually  inevitable.  Further,  if 
one  attempts  to  move  the  So  source  plane  farther  and  farther  away  from  the  data  plane  -  proposed 
here  for  isolating  physical  sources  -  problematical  ill-conditioning  is  to  be  expected.  Figure  1 1 
illustrates  a  simulation  case  in  which  spacings  that  were  not  at  all  extreme  produced  substantial 
noise  amplification  via  ill-conditioning. 

In  computational  practice  the  matrices  to  be  dealt  with  do  not  normally  present  themselves 
in  such  a  convenient  form  as  that  in  (21).  While  the  actual  matrices  from  (12)  are  full,  as  in  (22), 
their  insufficiency  is  much  less  evident.  They  are  much  larger;  the  full  matrix  system  is  not 
readily  seen  as  a  geometrical  rotation  of  a  clearer  system;  and  problematical  components  are  not 
readily  visible  as  being  too  similar  to  one  another.  At  the  same  time,  generalized  “rotations”  may 
be  applied  so  that  the  system  is  translated  into  an  eigenvalue-eigenvector  basis.  The  material 
below  pursues  ways  to  use  eigenvector  and  other  decompositions  of  the  system  to  restrict  the 
solution  to  well-supported  components. 

The  overall  strategy  here  of  evading  ill-conditioning  can  be  outlined  as  follows.  Consider  a 
linear  system 


L[q]  =  H 


(23) 
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in  which  L  is  some  linear  operator,  e.g.  the  matrix  implied  by  (12).  constitutes  the  data  and  we 
are  to  solve  for  q.  Re-express  the  quantities  H  and  q  so  that  they  are  expressed  by  series 
containing  a  new  basis. 


L 


j 


’LhjJv) 


(24) 


The  functions  and  ^  are  members  of  the  new  basis  and  and  q^  are  coefficients.  The  bases 
need  not  be  the  same  for  both  H  and  q  though  in  what  follows  they  will  always  be  so. 

Multiplying  (24)  by  any  one  basis  function  and  integrating  (taking  a  suitably  defined 
scalar  product  <...>)  produces  the  relation 


[«>,])  =  P5) 

j  m 

and  one  can  now  solve  for  the  coefficients  q^.  The  change  of  basis  is  selected  in  such  a  way  that 
the  system  can  be  truncated  sensibly,  altogether  deleting  equations  and  corresponding 
(transformed)  unknowns  that  have  weak  support.  The  basis  functions  are  constructed  so  that  they 
are  independent  and,  ideally,  orthogonal  under  the  scalar  product  definition,  i.e. 

(2«) 

where  is  the  Kronecker  delta.  Strong  independence  in  the  form  of  the  orthogonality  is 
advantageous  in  a  number  of  ways.  First,  we  noted  above  that  when  source  or  observation  points 
are  physically  close  to  one  another  so  that  their  contributions  to  the  system  are  similar,  ill- 
conditioning  will  tend  to  result.  Hence  here  we  seek  a  mathematical  space  in  which  to  operate  in 
which  the  known  and  unknown  quantities  are  associated  with  maximally  independent  influences 

in  the  system.  Secondly,  (26)  simplifies  computations  greatly.  In  matrix  terms, 

becomes  a  diagonal  matrix  and  therefore  no  onerous  manipulations  are  required  to  invert  the 
system.  Lastly  but  perhaps  most  important,  one  can  examine  the  elements  on  the  diagonal  matrix 
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directly  and  easily  which  of  them  are  weak,  hence  where  to  truncate  the  series  in 

(25). 


In  what  follows  we  examine  two  approaches  to  transforming  the  system  into  a  form  in 
which  weak  elements  are  discernible  and  readily  subject  to  elimination,  namely  use  of 
eigenvector  bases  and  of  sinusoidal  bases. 


5.2.1.  Eigenvector  Decomposition  and  Exclusion  of  Weak 
Eigenvalues 

The  eigenvector  decomposition  approach  to  dealing  with  ill-conditioning  has  been 
recognized  in  a  number  of  numerical  applications,  e.g.  [10].  One  proceeds  as  above,  using  (12)  to 
produce  an  algebraic  equation  of  the  form 

Aq  =  H  (27) 

where  7/  is  a  vector  of  the  data  and  the  vector  q  contains  the  point  sources  to  be  obtained.  The 
transformations  operate  directly  on  the  already  discretized  system,  not  on  the  original  continuous 
integral  relation  (11). 

The  matrix  A  is  decomposed  in  terms  of  its  eigenvectors  and  eigenvalues  as 

A  =  QAQ^  (28) 

where  g  contains  the  normalized  eigenvectors  and  A  is  a  diagonal  matrix  of  eigenvalues.  While 
similar  formulations  may  be  produced  for  rectangular  systems  of  equations,  in  all  examples  here 
square  matrices  are  used  so  we  remain  with  the  form  in  (28).  By  virtue  of  their  mutual 
orthogonality  property,  the  normalized  eigenvectors  form  a  matrix  with  the  property  that  its 
transverse  is  its  inverse: 

ee"  =  I  (29) 
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The  quantities  H  and  q  can  be  re-expressed  in  an  eigenvector  basis  as 


=  Qv 
=  Q 


Combining  these  equations  produces  an  algebraic  system  in  terms  of  quantities  in  the  transformed 
space,  i.e.  in  the  new  variables  77  =  H,  v  =  q  . 


7j  =  Av 


(31) 


The  inversion  of  this  system  is  immediate,  as  the  matrix  A  is  diagonal.  This  set  of  steps  is  the 
equivalent  of  starting  with  the  matrix  in  (22),  for  example,  and  transforming  it  to  that  in  (21). 
Ultimately,  source  values  q  can  be  obtained  from  the  solution  for  the  transformed  variable  v  via 
(30). 


The  point  of  these  manipulations  is  that,  as  in  our  treatment  of  (21),  one  can  discern  readily 
in  the  transformed  system  what  equations  are  too  weak  to  merit  inclusion  in  the  total  system. 
Equations  in  (31)  containing  eigenvalues  that  are  smaller  than  some  cutoff  a  relative  to  the 
maximum  (absolute  value)  X-max  are  discarded.  That  is,  the  complete  solution  v  is  approximated 
by  V, 


y.  AAmax  >  « 
<  A- 


(32) 


from  which 


q  ~ 


(33) 


Note  that  this  produces  the  same  number  of  elements  in  q  as  before  (72),  but  simply  using  a  subset 
of  m  of  the  complete  set  of  n  eigenvectors.  Using  the  complete  set  of  eigenvectors  {m  =  n)  will 
produce  strict  equality  between  q  and 
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To  investigate  and  to  illustrate  eigenvector  decomposition  we  pursue  some  simple  1-D 
examples.  The  solution  is  taken  to  be  ^  =  sin  /  33)  and  a  100  x  100  symmetric  matrix  A  is 

produced  using  a  random  number  generator,  resulting  in  C  ~  5x10^.  H  is  obtained  from^^. 
Error  in  the  data  H  is  produced  by  adding  random  white  noise  over  a  band  of  values  \-b,  b\,  where 


RelNoiseMag 


2b 


0.05 


(34) 


This  combination  of  C  and  RelNoiseMag  will  normally  produce  an  unacceptably  distorted 
solution:  the  SNR  (loosely  defined)  is  on  the  order  of  20,  while  something  on  the  order  of  5000  or 
better  is  needed.  Figure  26  shows  the  right  hand  side  (data,  H)  obtained  from  the  A  and  q  that 
have  been  specified.  Note  that  the  noise  is  considerably  less  than  is  typically  the  case  in  real  field 
data. 


Figure  26.  Right  hand  side  H  in  (27),  with  and  without  noise  at  RelNoiseMag  =  5%. 

Without  added  noise,  the  precision  of  the  data  H  (roundoff  error  in  the  computer)  still 
securely  places  its  SNR  above  that  required  by  the  condition  number.  Hence  the  noise  free 
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solution  of  (27),  not  shown,  is  extremely  accurate,  with  no  graphically  visible  distortions.  At  the 
other  extreme,  direct  solution  of  the  system  with  added  noise  at  the  indicated  level  produces  q 
values  that  are  predominately  just  amplified  noise.  If  one  retains  the  noise-free  data  and  truncates 
the  eigen-system,  using  a  cutoff  of  a  =  10'^  with  resulting  min  =  85/100,  the  solution  in  Figure 
27,  top,  results.  The  solution  visibly  follows  the  actual  (exact)  solution.  However,  truncating  the 
eigen-system  at  a  “safe”  level  discards  enough  meaningful  information  so  that,  in  translation  back 
from  vto  q,  the  deletions  cause  lumps  and  bumps.  This  highlights  a  basic  problem  in  the  eigen- 
system  truncation  approach.  Namely,  the  nature  of  different  portions  of  the  information  in  the 
eigen-system  may  not  mean  that  retaining  values  above  a  certain  level  of  significance  will 
translate  into  values  above  that  level  of  significance  in  the  original  space. 

Figure  27,  bottom,  shows  the  solution  obtained  by  the  same  eigen-system  truncation  but 
using  the  noisy  data  in  Figure  26.  There  are  discernible  differences  between  the  solution  from  the 
noise-free  and  noisy  data.  However  these  are  slight.  The  eigen-system  truncation  has  indeed 
restrained  the  influence  of  noise  in  the  data  to  a  minor  factor,  although  the  truncation  itself  has 
still  left  some  distortions. 
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Figure  27.  Top:  truncated  eigen-solution  without  noise  added  to  data.  Bottom:  same,  but 
when  noise  was  added  to  the  data. 
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5.2.2.  Fourier  Decomposition  and  Exciusion  of 
Insufficiency  Supported  Spatiai  Frequencies 

The  formulation  in  this  approach  operates  directly  on  the  continuous  variables  and  integral 
expression  (11),  producing  an  alternative  discrete  system  which  is  readily  truncated  to  advantage. 
It  exploits  relations  and  strategies  outlined  in  [2].  Essentially,  the  system  is  translated  into  the 
(spatial)  frequency  domain,  in  terms  of  spatial  2-D  Fourier  transforms  of  both  the  data  and  of  the 
source  distribution.  This  allows  one  to  solve  for  coefficients,  again  without  in  fact  inverting 
any  matrices,  truncating  the  system  at  spatial  wave  numbers  beyond  which  contributions  are 
negligible.  Note  that,  without  treatment,  the  solution  noise  that  “slips  through”  the  linear  operator 
from  the  data  tends  to  correspond  to  high  spatial  frequencies.  In  Figure  1 1,  the  very  substantial 
distortions  in  the  solution  appear  in  the  form  of  values  oscillating  between  relatively  large 
positive  and  negative  values  between  neighboring  points.  From  one  point  of  view,  this  is  a 
relatively  benign  form  of  noise  in  that  it  tends  to  integrate  out  when  one  uses  (1 1)  to  calculate 
field  values  at  some  elevation.  The  (numerically  expressed)  integral  has  something  of  a  “blind 
spot”  for  information  at  this  frequency.  This  is  convenient  for  upward  continuation.  However, 
when  one  proceeds  in  the  inverse  direction,  the  integral  and  its  numerical  equivalent  are  not  good 
at  restricting  data  noise  of  this  wavelength  from  assuming  an  arbitrary  large  value  in  the  solution. 
Our  strategy  will  be  to  recast  the  system  in  the  frequency  domain  and  then  to  forcibly  truncate  the 
participation  of  frequencies  higher  than  the  system  will  readily  support  accurately. 

To  proceed,  express  as 

M  N 

«'=(>■’)  =  (35) 

m=l  n=l 


Equation  (12)  becomes 


(>.(>•)  =  Z /‘(•s 


^  B.  ■(  (x- x') i  +  (y-  y )y  +  o' z) 


[x-x')^  +{y-y')^ 


n3/2 


(36) 
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where,  conveniently,  we  benefit  from  the  fact  that  (z  -  z')  is  uniformly  equal  to  some  depth  d 

(both  data  and  sources  are  assumed  to  be  located  on  planes  of  constant  z).  In  principle,  the 
expression  in  (35)  could  be  written  as  an  integral  of  continuous  quantities  rather  than  a  finite  sum 
of  discrete  terms.  In  that  case,  any  two  different  k  values,  say,  ki  and  k^  will  be  associated  with 
basis  functions  that  fulfill  our  orthogonality  requirement  that 


—00 


(37) 


SO  that  (26)  applies.  In  numerical  practice,  the  integration  is  over  some  suitable  finite  interval  and 
a  corresponding  finite  selection  of  k  values  is  used,  as  in  the  summation  (35).  For  example, 
driven  either  by  the  need  for  a  certain  minimum  domain  size  or  a  certain  spacing  of 
wavenumbers,  one  might  pick  the  domain  size  L  with  k^  =  InmlL.  Matlab,  in  which  language 
system  the  computations  here  were  performed,  automatically  applies  a  k  selection  of  this  sort  in 
its  FFTs  and  IFFTs  so  that  a  discrete  form  of  (37)  holds  exactly. 


With  a  change  of  variables 

u  =  x-x\  v  =  y-y'  (38) 

(36)  becomes 


ou 

*,(>•)  =  Z  J  dS 
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{u’f+{v'f+d^ 


3/2 


(39) 


where  V  need  only  be  calculated  once  for  each  d.  It  is  basically  just  a  weighted  FT  of  the  Green 
function  for  referred  to  the  origin.  A  significant  attraction  of  this  approach  is  that  all  r 

dependence  of  (beyond  d)  is  contained  in  the  exponentials  outside  the  integral. 

Values  of  V  for  each  d  can  simply  be  stored  and  used  for  arbitrary  r  as  required. 
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The  remaining  exponentials  in  (39)  can  be  integrated  out  if  one  expresses  itself  as 


m',n' 

for  the  same  set  of  wavenumbers  as  were  used  for  q^.  The  coefficients  p^,^,  are  simply  obtained 

from  the  FT  of  the  data.  To  conclude,  substitute  (40)  into  (39),  multiply  by  and 

integrate  (x,y)  each  over  (^— oo,qo)  or  over  the  aforementioned  suitable  interval,  to  obtain 

Aj  =  V(k,,k.)  (41) 

for  each  {ij).  As  in  the  previous  approach,  all  quantities  are  now  scalars  and  no  matrix  solution  is 
required. 

The  truncation  is  accomplished  by  retaining  only  spatial  frequencies  that  contribute  to  the 
solution  securely: 


>  a 

yiy.^  s  a  ^ 


(42) 


While  geographically  applicable  may  make  the  k  dependence  of  V  non-monotonic,  the  system 

in  (42)  succeeds  in  filtering  out  all  frequencies  with  insufficient  support,  whether  high  or  low.  In 
any  case,  the  nature  of  Fourier  decomposition  is  such  that  k  values  above  some  level  will 
inevitably  contribute  little  to  both  the  data  and  the  solution,  whatever  the  pattern  of  dependence 
on  lower  k  values.  High  frequency  noise  influences  will  be  filtered  out. 

Equation  (41)  indicates  that  individual  Fourier  components  of  the  data  are  only  associated 
with  individual,  corresponding  Fourier  components  of  the  source  distribution.  This  prompts  the 
question  of  whether  the  necessary  truncation  of  k  space  for  V  will  be  sufficient  to  define  well. 
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This  relates  in  part  to  the  frequency  resolution  level  that  is  inherent  in  the  integral  operator. 
Probably  of  greater  concern  is  the  specific  frequency  content  of  the  data,  i.e.  whether  the  range  of 
significantly  contributing  wavenumbers  in  the  data  is  sufficient  for  the  requirements  of  aij.  The 
spatial  distribution  of  the  data  Bq  will  in  general  be  broader  than  that  of  as  deeper  and  deeper  Sq 
are  considered  and  becomes  increasingly  concentrated.  Higher  and  higher  k  will  be  required 
for  which  may  not  be  supported  by  Bq,  i.e.  the  corresponding  j  may  be  essentially  noise. 

To  explore  these  questions  and  to  test  the  frequency  truncation  approach  overall,  we  begin 
with  a  simplified  simulation  case  in  which  the  earth’s  field  ( )  is  vertical  and  the  responding 

subsurface  body  is  a  vertical  dipole.  This  produces  symmetrical  signal  patterns  that  are  easy  to 
assess  visually.  With  the  target  2itd=  -0.5  m  and  a  =  0.01,  the  solutions  in  Figure  28  result.  The 
upper  right  plot  in  the  figure  shows  what  perfect,  noise-free  data  would  look  like.  The  two  plots 
of  the  noisy  data  reconstructed  from  the  truncated  q^  solutions  (middle  and  bottom  right) 
illustrate  the  level  of  clutter  in  the  noisy  data.  Without  either  eigen-system  or  frequency  value 
truncation,  a  direct  solution  for  q  contains  very  large  oscillations  between  neighboring  points 
(upper  left  plot).  In  contrast,  both  of  the  truncated  systems  produce  tolerable  renderings  of  an 
equivalent  source  distribution  half  a  meter  above  the  actual,  infinitely  concentrated  physical 
source  (middle  and  bottom  left  plots). 

Lowering  the  source  place  Sq  in  order  to  cause  the  q  pattern  to  focus  on  the  physical  source 
location  produces  the  results  in  Figure  29.  The  truncated  solutions  succeed  at  least  in  the  sense 
that  they  remain  coherent,  with  noise  suppressed.  By  contrast,  a  direct  q  solution  (no  truncation, 
not  shown)  for  Zq  =  -0.5  m  covers  whole  plane,  oscillating  between  +  10^^.  However  neither  of 
the  truncated  solutions  focuses  notably  on  the  infinitesimal  physical  source.  As  Zq  drops  below 
zero  (“ground  surface”)  and  the  condition  number  rises  exponentially,  the  eigen-system 
truncation  allows  fewer  and  fewer  eigenvectors  to  participate  {min  decreases,  see  plot  headers). 
Similarly,  the  resolution  of  the  frequency  truncation  system  is  also  limited,  either  because  of  the 
limited  range  of  the  spectrum  of  the  data  or  of  the  (numerical)  integral  operator,  or  both.  The 
finite  focusing  of  the  truncated  solutions  in  the  figures  is  summed  up  in  Table  1.  These  results 
suggest  that,  while  both  methods  stabilize  the  solution  impressively  in  the  face  of  enormous 
condition  numbers,  neither  is  likely  to  be  useful  for  focusing  on  subsurface  source  locations. 
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Figure  28.  For  data  taken  20  cm  above  the  surface,  with  the  source  plane  5^0  on  the  ground 
surface  (Zq  =  0  cm).  Left:  Direct  q  solution  (top);  that  from  the  truncated  FT 
(middle);  and  that  from  the  truncated  eigen-system  (bottom).  Right:  Data  without 
added  noise  (top);  noisy  data  reconstructed  from  truncated  FT  source  solution 
(middle);  same,  reconstructed  from  truncated  eigen-system  (bottom). 
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Figure  29.  Source  solutions  for  Sq  surfaces  deeper  than  that  (Zq  =  0)  in  Figure  28.  Here  Zq 
=  -0.25  (top)  and  -0.5  (bottom),  the  latter  being  the  depth  of  the  point  physical 
source.  Left:  q  from  frequency  truncation.  Right:  from  eigen-system 
truncation. 


Table  1.  Depth  Zq  of  surface  So  containing  equivalent  sources  q;  approximate  spatial  width 
of  the  truncated  q  solution;  condition  number  C  of  the  governing  matrix  before 
truncation. 


Zq 

approx  q  width 

C 

0 

1  m 

~  10^ 

-0.25  m 

0.5  m 

~4x  10*^ 

-0.5  m 

0.6  m 

~2x  10^^ 
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To  distinguish  better  between  the  performances  of  the  two  truncation  systems,  we  perform 
another  test  with  somewhat  more  realistic  i.e.  complex  patterns  in  the  relevant  quantities.  While 
the  target  is  again  at  Z  =  -0.5  m,  its  moment  here  is  m  =  [0,1,1],  as  if  it  is  tilted  (the  previous 

vertical  dipole  had  m  =  [0,0,1] ).  Further,  a  non- vertical  B^is  assumed,  in  particular,  that  which 

applied  at  the  APG  test  site,  =  [-0.08,  0.38,  -0.92].  Computations  for  both  truncation  systems 

were  performed  with  a  =  0.01  and  RelNoiseMag  =  0.01.  The  changes  introduced  relative  to  the 
first  test  produce  more  complicated  signal  distributions  and  F(k)  patterns.  Evidently  the  eigen- 
structure  of  the  system  is  also  more  complicated  in  such  a  way  that  simple  truncation  has 
unanticipated  and  unwelcome  effects.  The  computed  q  distributions  in  Figure  30  show  that 
frequency  truncation  produced  an  orderly  shape  reflecting  the  Y-Z  tilt  of  the  target.  The  eigen- 
system  truncation  results  are  irregular  and  do  not  show  the  physical  Y-Z  spread  of  source  action 
in  the  central  cluster.  Further  (Figure  31),  when  the  data  on  the  measurement  plane  is 
reconstructed  using  these  truncated  q  solutions,  that  from  the  frequency  truncation  approach  is 
much  superior  in  both  magnitude  and  shape. 

Whether  this  picture  might  be  altered  by  some  sort  of  strategic  manipulation  of  particulars 
in  the  truncation  systems  -  e.g.  adjustment  of  a  or  other  means  of  selecting  residual  participants 
in  the  eigen-solution  -  must  be  determined  by  subsequent  research.  Overall,  as  they  stand  these 
results  suggest  that  eigen-system  truncation  could  not  be  relied  upon  either  for  downward 
solution  for  source  concentrations  or  for  computation  of  above-ground  signals.  Frequency 
truncation  appears  able  to  produce  tolerable  equivalent  source  distributions,  albeit  limited  in  their 
concentration,  with  the  ability  to  compute  consistent  above-ground  signal  patterns. 
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Figure  30.  Truncated  solutions  for  equivalent  source  distributions  with  Zq  =  -0.25  when  the 
physical  target  is  at  -0.5  m.  Top:  frequency  truncation.  Bottom:  Eigen-system 
truncation. 
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Figure  31.  Distributions  of  the  data  on  the  measurement  plane  Top:  “Actual”  data  with 
(left)  and  without  noise  (right).  Bottom:  Data  on  measurement  plane 
reconstructed  from  truncation  solutions  based  on  the  noisy  data,  derived  from  the 
frequency  (left)  and  eigenvalue  truncation  solutions  for  q  (right). 
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5.2.3.  Control  of  Clutter  Amplification  as  Well? 


As  a  last  experiment,  consider  a  case  again  with  the  earth’s  field  vertical,  including  a 
vertical  dipole  target  (  m  =  [0,  0,  1] )  at  60  cm  depth,  offset  from  which  is  a  piece  of  clutter  (  m  = 
[0,  0,  10'^]  at  20  cm  depth.  In  order  to  focus  on  amplification  of  the  clutter  source  at  depth  (see 
Section  ***),  we  assume  no  random  noise  in  the  data  such  might  be  amplified  by  any  ill- 
conditioning.  As  So  is  lowered  the  computations  produce  the  equivalent  source  distributions  in 
Figure  32,  both  with  (a  =  0.01)  and  without  frequency  truncation.  For  a  source  plane  Sq  on  the 
“ground  surface”  (Zq  =  0)  above  both  items,  both  direct  (non-truncated)  and  truncated  solutions 
are  the  same  .  At  Zq  equal  to  the  clutter  depth,  the  solutions  are  dominated  appropriately  by  the 
clutter  item,  more  so  in  the  direct  solution.  With  deeper  Sq  the  direct  solution  is  increasingly 
dominated  by  the  sources  pertaining  to  the  clutter,  which  spread  and  increase  in  magnitude. 
Additionally,  with  C  increasing  to  ~  5  x  10^^  at  Zq  =  60  cm,  any  noise  whatsoever  in  the  data 
would  have  further  swamped  the  direct  solution,  causing  high  amplitude,  high  frequency 
oscillations  over  the  entire  domain.  However,  in  our  noiseless  example,  the  frequency-truncated 
solution  retains  a  more  appropriate  character  and  magnitude,  some  distortion  notwithstanding. 
Perhaps  most  notably,  the  frequency  truncation  appears  to  restrain  the  tendency  of  the  clutter- 
related  sources  to  dominate  the  solution  at  depths  greater  than  the  clutter  item  itself  This  implies 
some  limited  deletion  of  signal  content  in  the  sense  that  a  reconstruction  of  the  data  on  the 
measurement  plane  from  the  truncated  ^  at  Zq  =  60  cm  shows  limited  resolution  (Figure  33). 
However,  the  gross  features  and  appropriate  magnitude  of  the  data  are  retained. 
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Figure  32.  Equivalent  source  distributions  computed  on  planes  at  various  depths,  Zq.  Left: 
Direct  solution  (no  truncation).  Right:  Frequency  truncation. 
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Figure  33..  Same  case  as  above,  showing  reconstruction  of  the  data  on  the  measurement 
plane  from  the  truncated  g  at  Zq  =  60  cm. 
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